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Abstract 

In  contrast  to  discrete  descriptions  of  fracture,  phase-field  descriptions  do  not  require  numerical  tracking 
of  discontinuities  in  the  displacement  field.  This  greatly  reduces  implementation  complexity.  In  this  work, 
we  extend  a  phase-field  model  for  quasi- static  brittle  fracture  to  the  dynamic  case.  We  introduce  a  phase- 
field  approximation  to  the  Lagrangian  for  discrete  fracture  problems  and  derive  the  coupled  system  of 
equations  that  govern  the  motion  of  the  body  and  evolution  of  the  phase-field.  We  study  the  behavior  of 
the  model  in  one  dimension  and  show  how  it  influences  material  properties.  For  the  temporal  discretization 
of  the  equations  of  motion,  we  present  both  a  monolithic  and  staggered  time  integration  scheme.  We 
study  the  behavior  of  the  dynamic  model  by  performing  a  number  of  two  and  three  dimensional  numerical 
experiments.  We  also  introduce  a  local  adaptive  refinement  strategy  and  study  its  performance  in  the  context 
of  locally  refined  T- splines.  We  show  that  the  combination  of  the  phase-field  model  and  local  adaptive 
refinement  provides  an  effective  method  for  simulating  fracture  in  three  dimensions. 

Keywords:  Phase  field,  Fracture  mechanics,  Isogeometric  analysis,  Adaptive  refinement,  T-splines 


1  Introduction 

The  prevention  of  fracture-induced  failure  is  a  major  constraint  in  engineering  designs,  and  the  numerical 
simulation  of  fracture  processes  often  plays  a  key  role  in  design  decisions.  As  a  consequence,  a  wide  variety 
of  fracture  models  have  been  proposed.  A  particularly  successful  model  is  provided  by  Griffith’s  theory  for 
brittle  fracture,  which  relates  crack  nucleation  and  propagation  to  a  critical  value  of  the  energy  release  rate. 
A  general  concept  in  the  Griffith’s -type  brittle  fracture  models  is  that  upon  the  violation  of  the  critical  energy 
release  rate  a  fully  opened  crack  is  nucleated  or  propagated.  As  a  consequence,  the  process  zone,  i.e.,  the 

*  Submitted  to  Computer  Methods  in  Applied  Mechanics  and  Engineering. 
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zone  in  which  the  material  transitions  from  the  undamaged  to  the  damaged  state  is  lumped  into  a  single  point 
at  the  crack  tip. 

Due  to  the  complexity  of  fracture  processes  in  engineering  applications,  numerical  methods  play  a  crucial 
role  in  fracture  analyses.  In  particular,  finite  element  methods  are  used  extensively  in  conjunction  with 
Griffith’ s-type  linear  elastic  fracture  mechanics  models.  Among  the  most  commonly  used  finite  element 
models  are  the  virtual  crack  closure  technique  (see  Krueger  (2004)),  and,  in  more  recent  years,  the  extended 
finite  element  method  introduced  by  Moes,  Dolbow,  and  Belytschko  (1999).  All  of  these  approaches  represent 
cracks  as  discrete  discontinuities,  either  by  inserting  discontinuity  lines  by  means  of  remeshing  strategies,  or 
by  enriching  the  displacement  field  with  discontinuities  using  the  partition  of  unity  method  of  Babuska  and 
Melenk  (1997).  Tracing  the  evolution  of  complex  fracture  surfaces  has,  however,  proven  to  be  a  tedious  task, 
particularly  in  three  dimensions. 

Recently,  alternative  methods  for  the  numerical  simulation  of  brittle  fracture  have  emerged.  In  these 
approaches,  discontinuities  are  not  introduced  into  the  solid.  Instead,  the  fracture  surface  is  approximated  by 
a  phase-field,  which  smoothes  the  boundary  of  the  crack  over  a  small  region.  The  major  advantage  of  using  a 
phase-field  is  that  the  evolution  of  fracture  surfaces  follows  from  the  solution  of  a  coupled  system  of  partial 
differential  equations.  Implementation  does  not  require  the  fracture  surfaces  to  be  tracked  algorithmically. 
This  is  in  contrast  to  the  complexity  of  many  discrete  fracture  models,  and  is  anticipated  to  be  particularly 
advantageous  when  multiple  branching  and  merging  cracks  are  considered  in  three  dimensions. 

A  phase-field  model  for  quasi-static  brittle  fracture  emanated  from  the  work  of  Bourdin  and  co-workers 
on  the  variational  formulation  for  Griffith’ s-type  fracture  models  (see  Bourdin,  Francfort,  and  Marigo  (2008) 
for  a  comprehensive  overview).  The  variational  formulation  for  quasi- static  brittle  fracture  leads  to  an  energy 
functional  that  closely  resembles  the  potential  presented  by  Mumford  and  Shah  (1989),  which  is  encountered 
in  image  segmentation.  A  phase-field  approximation  of  the  Mumford-Shah  potential,  based  on  the  theory 
of  T-convergence,  was  presented  by  Ambrosio  and  Tortorelli  (1990).  This  approximation  was  adopted  by 
Bourdin,  Francfort,  and  Marigo  (2008)  to  facilitate  the  numerical  solution  of  their  variational  formulation. 
Recently,  this  model  has  been  applied  in  a  dynamic  setting  by  Bourdin,  Larsen,  and  Richardson  (2011); 
Larsen,  Ortner,  and  Siili  (2010);  and  Larsen  (2010),  but  application  to  structures  of  engineering  interest  has 
not  been  considered. 

An  alternative  quasi-static  formulation  of  this  phase-field  approximation  has  been  presented  in  the  re¬ 
cent  work  of  Miehe,  Hofacker,  and  Welschinger  (2010a)  and  Miehe,  Welschinger,  and  Hofacker  (2010b). 
In  this  formulation,  the  phase-field  approximation  follows  from  continuum  mechanics  and  thermodynamic 
arguments.  Besides  the  provision  of  an  alternative  derivation,  Miehe,  Hofacker,  and  Welschinger  (2010a) 
also  added  various  features  to  the  model  that  are  key  to  its  application  to  engineering  structures. 

Independently  from  the  phase-field  formulation  based  on  Griffith’s  theory,  dynamic  phase-field  fracture 
models  have  been  developed  based  on  Landau-Ginzburg  type  phase-field  evolution  equations,  e.g.,  Karma, 
Kessler,  and  Levine  (2001).  However,  we  favor  the  phase-field  formulation  of  the  Bourdin-type  since  the 
physical  properties  of  Griffith’s  theory  are  well  understood  and  have  proven  useful  in  engineering  applica¬ 
tions. 

In  this  contribution  we  extend  the  quasi-static  model  presented  by  Miehe,  Hofacker,  and  Welschinger 
(2010a)  to  the  dynamic  case.  We  begin  by  formulating  the  Lagrangian  for  the  discrete  fracture  problem 
in  terms  of  the  displacements  and  the  phase-field  approximation  of  the  crack  path.  Then,  using  the  Euler- 
Lagrange  equations  of  this  approximation,  we  derive  the  strong  form  equations  of  motion.  A  detailed  analysis 
of  the  analytical  one-dimensional  solution  to  the  strong  form  equations  is  then  presented.  In  this  analysis  we 
show  that  although  the  length- scale  parameter  associated  with  the  phase-field  approximation  is  introduced 
as  a  numerical  parameter  it  is,  in  fact,  a  material  parameter  that  influences  the  critical  stress  at  which  crack 
nucleation  occurs. 

The  numerical  solution  of  the  strong  form  of  the  equations  of  motion  requires  a  spatial  and  temporal  dis¬ 
cretization.  We  formulate  the  spatial  discretization  by  means  of  the  Galerkin  method.  In  this  work,  we  have 
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used  NURBS  and  T-spline  basis  functions  as  the  finite  dimensional  approximations  to  the  function  spaces 
of  the  weak  form.  However,  we  note  that  standard  C° -continuous  finite  elements  could  also  be  used.  For 
the  temporal  discretization,  we  present  monolithic  and  staggered  time  integration  schemes.  The  monolithic 
scheme  requires  solution  of  the  coupled  equations  simultaneously.  For  the  staggered  scheme,  the  displace¬ 
ments  (via  the  momentum  equation)  and  phase-field  (via  the  phase-field  equation)  are  solved  for  separately. 
This  provides  flexibility  in  solution  strategies  by  allowing  the  momentum  equations  to  be  solved  either  im¬ 
plicitly  or  explicitly. 

To  conclude  this  paper,  we  study  the  behavior  of  the  model  by  performing  a  number  of  numerical  bench¬ 
mark  experiments  for  crack  propagation  and  branching.  These  experiments  show  that  the  phase-field  model 
can  capture  complex  crack  behavior  in  both  two  and  three  dimensions,  without  introducing  any  ad  hoc  crite¬ 
ria  for  crack  nucleation  and  branching.  In  addition,  we  propose  an  adaptive  refinement  strategy  that  allows 
for  the  efficient  simulation  of  complex  crack  patterns.  We  show  that  adaptive  refinement  maintains  accuracy 
while  providing  greater  efficiency  in  terms  of  the  number  of  degrees-of-freedom.  Finally,  we  apply  the  adap¬ 
tive  refinement  strategy  to  a  three-dimensional  problem.  This  final  problem  illustrates  the  potential  strength 
of  the  phase-field  model,  which  is  the  ability  to  efficiently  model  dynamic  fracture  in  three-dimensions. 


2  Formulation 

We  briefly  introduce  Griffith’s  theory  for  dynamic  brittle  fracture  in  bodies  with  arbitrarily  discrete  cracks. 
We  then  present  a  phase-field  formulation  as  a  continuous  approximation  of  the  discrete  fracture  model.  We 
conclude  this  section  with  a  study  of  analytical  solutions  of  the  phase-field  in  the  one-dimensional  quasi-static 
case,  which  reveals  many  interesting  features  of  the  model. 


2.1  Griffith’s  theory  of  brittle  fracture 


We  consider  an  arbitrary  body  U  C  (with  d  G  {1,2,  3})  with  external  boundary  dfl  and  internal  dis¬ 
continuity  boundary  T  (see  Figure  la).  The  displacement  of  a  point  x  G  U  at  time  t  G  [0,  T]  is  denoted  by 
u(x,  t)  G  Md.  Spatial  components  of  vectors  and  tensors  are  indexed  by  i,  j  =  1, . . . ,  d.  The  displacement 
field  satisfies  time-dependent  Dirichlet  boundary  conditions,  Ui(x ,  t)  =  gi(x ,  t),  on  dClg.  C  cftl,  and  time- 
dependent  Neumann  boundary  conditions  on  dGthi  C  d£l.  We  assume  small  deformations  and  deformation 
gradients,  and  define  the  infinitesimal  strain  tensor,  e(x,t)  G  Rdxd ,  with  components 


_  1 

£ij  ~  U(iJ)  ~  2 


f  duj  \ 

\  dxj  dxi  J 


(1) 


as  an  appropriate  deformation  measure.  We  assume  isotropic  linear  elasticity,  such  that  the  elastic  energy 
density  is  given  by 

=  (2) 

with  A  and  g  the  Lame  constants.  Note  that  we  use  the  Einstein  summation  convention  on  repeated  indices. 

The  evolving  internal  discontinuity  boundary,  T(t),  represents  a  set  of  discrete  cracks.  In  accordance  with 
Griffith’s  theory  of  brittle  fracture,  the  energy  required  to  create  a  unit  area  of  fracture  surface  is  equal  to  the 
critical  fracture  energy  density  Qcl .  The  total  potential  energy  of  the  body,  ^ pot ,  being  the  sum  of  the  elastic 
energy  and  the  fracture  energy,  is  then  given  by 


^ pot 


n  r 


Gcdx 


(3) 


1  This  critical  fracture  energy  density  is  commonly  referred  to  as  the  critical  energy  release  rate,  or,  in  the  context  of  cohesive  zone 
models,  the  fracture  toughness. 
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Figure  1:  (a)  Schematic  representation  of  a  solid  body  Cl  with  internal  discontinuity  boundaries  T.  (b)  Ap¬ 
proximation  of  the  internal  discontinuity  boundaries  by  the  phase-field  c(x,  t).  The  model  parameter  e  con¬ 
trols  the  width  of  the  failure  zone. 


where  we  have  defined  the  symmetric  gradient  operator,  Vs  :  u  e,  as  a  mapping  from  the  displacement 
field  to  the  strain  field.  Since  brittle  fracture  is  assumed,  the  fracture  energy  contribution  is  merely  the  critical 
fracture  energy  density  integrated  over  the  fracture  surface.  In  the  case  of  small  deformations,  irreversibility 
of  the  fracture  process  dictates  that  T(t)  C  T(t  +  At)  for  all  At  >  0.  Hence,  translation  of  cracks  through 
the  domain  is  prohibited,  but  cracks  can  extend,  branch,  and  merge. 

The  kinetic  energy  of  the  body  Q  is  given  by 


^ kin 


pUiliidx 


(4) 


with  um  ^  and  p  the  mass  density  of  the  material.  Combined  with  the  potential  energy  (3)  this  renders  the 
Lagrangian  for  the  discrete  fracture  problem  as 


L(u,U,T)  =  ^kin(u)  -^pot(u,T)  =  J 


-piiiUi  ~i)e(ysu) 


d x  —  J  Gc^x 

r 


(5) 


The  Euler-Lagrange  equations  of  this  functional  determine  the  motion  of  the  body.  From  a  numerical  stand¬ 
point,  tracking  the  evolving  discontinuity  boundary,  T,  often  requires  complex  and  costly  computations.  This 
is  particularly  so  when  interactions  between  multiple  cracks  (even  in  two  dimensions),  or  complex  shaped 
cracks  in  three  dimensions  are  considered.  Of  particular  interest  in  the  case  of  dynamic  fracture  simulations, 
as  considered  in  this  work,  is  the  ability  to  robustly  model  crack  branching.  In  the  remainder  of  this  work  we 
pursue  a  formulation  which  is  capable  of  handling  cracks  of  arbitrary  topological  complexity. 


2.2  Phase-field  approximation 

In  order  to  circumvent  the  problems  associated  with  numerically  tracking  the  propagating  discontinuity  rep¬ 
resenting  a  crack,  we  approximate  the  fracture  surface,  T,  by  a  phase-field,  c(x,t)  E  [0, 1].  The  value  of 
this  phase-field  is  equal  to  1  away  from  the  crack  and  is  equal  to  0  inside  the  crack  (see  Figure  lb).  We 
employ  the  approximation  as  discussed  in  Miehe,  Hofacker,  and  Welschinger  (2010a),  which  is  essentially 
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an  extension  of  the  phase-field  approximation  introduced  by  Bourdin,  Francfort,  and  Marigo  (2008)2.  As  in 
Bourdin,  Francfort,  and  Marigo  (2008),  we  approximate  the  fracture  energy  by 


J  Gr  dx 

r 


/' 


(c— l)2  dc  dc 

4e  dxj  dxj 


dx 


(6) 


where  e  E  M+  is  a  model  parameter  that  controls  the  width  of  the  smooth  approximation  of  the  crack.  From 
equation  (6)  it  is  clear  that  a  crack  is  represented  by  regions  where  the  phase-field,  c(sc,  £),  goes  to  zero.  As 
elaborated  by  e.g.  Bourdin,  Francfort,  and  Marigo  (2008),  in  the  limit  of  the  length  scale  e  going  to  zero,  the 
phase-field  approximation  converges  to  the  discrete  fracture  surface. 

To  model  the  loss  of  material  stiffness  in  the  failure  zone  (i.e.,  the  phase-field  approximation  of  the 
fracture  surface),  the  elastic  energy  is  approximated  by 

ipe(e,c)  «  [(1  -  k)c2  +  A#+(e)  +^~(e).  (7) 


As  in  Miehe,  Hofacker,  and  Welschinger  (2010a),  we  distinguish  between  the  cases  of  compressive  and  ten¬ 
sile  loading.  By  only  applying  the  phase-field  parameter  to  the  tensile  part  of  the  elastic  energy  density, 
we  prohibit  crack  propagation  under  compression.  This  model  feature  has  been  observed  to  be  particularly 
important  in  dynamic  simulations,  as  stress  waves  reflecting  from  domain  boundaries  tend  to  create  phys¬ 
ically  unrealistic  fracture  patterns.  The  model  parameter  k  <C  1  is  introduced  to  prevent  the  positive  part 
of  the  elastic  energy  density  from  disappearing  when  the  phase-field  is  equal  to  zero,  which  has  been  ob¬ 
served  to  improve  computational  robustness  in  the  quasi-static  simulations  presented  by  Miehe,  Hofacker, 
and  Welschinger  (2010a). 

Substitution  of  the  phase-field  approximations  for  the  fracture  energy  (6)  and  the  elastic  energy  density 
(7)  into  the  Lagrange  energy  functional  (5)  yields 


Le(u,  ii,  c)  =  J  (^puiUi  -  [(1  -  k)c2  +  k]ip+(ysu)  -  ipe  (Vsu)^  d® 

(c  —  l)2  dc  dc 


-  f  Sc 

Jn 


4e 


+  e 


dxi  dxi 


(8) 


dx. 


Note  that  in  order  to  conserve  mass  the  kinetic  energy  term  is  unaffected  by  the  phase-field  approximation. 
The  dependence  of  the  Lagrange  energy  functional  on  the  propagating  discontinuity  boundary  is  now  captured 
by  the  phase-field,  c(x,t)9  which  simplifies  the  numerical  treatment  of  the  model.  In  Miehe,  Hofacker, 
and  Welschinger  (2010a)  an  additional  viscosity  contribution  is  introduced.  For  the  dynamic  simulations 
performed  within  this  study  no  beneficial  effects  of  this  viscosity  parameter  were  encountered,  and  hence  this 
term  is  omitted  for  brevity. 

Now  that  we  have  formulated  the  Lagrangian  in  terms  of  the  independent  fields  u(x,  t )  and  c(x1  £),  we 
can  use  the  Euler-Lagrange  equations  to  arrive  at  the  strong  form  equations  of  motion 


< 


( 4e(l  -  fc)^>+ 

V  Sc 


da. 


*3 


c  —  4eJ 


dx 

d2 


=  pui 


dx2 


=  1 


on  fix]0,T[ 
on  fix]0,T[ 


(9) 


2For  the  most  part,  we  adopt  the  notation  introduced  by  Bourdin,  Francfort,  and  Marigo  (2008),  the  exception  being  the  use  of  k  in 
place  of  rj  for  the  model  parameter  that  controls  conditioning  of  the  linear  system.  In  Miehe,  Hofacker,  and  Welschinger  (2010a),  use  is 
made  of  a  damage  field  d(x,  t)  =  1  —  c(x,  t),  and  a  length  scale  l  =  2e  to  control  the  width  of  the  failure  zone. 
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where  u  =  4^  and  the  Cauchy  stress  tensor  a  £ 


ixd 


is  defined  by 


dipt  , 


=  [(i-^  +  ^  + 


(10) 


These  equations  of  motion  can  be  solved  to  find  both  the  displacement  field  u(x,t)  and  phase-field  c(x,  t ). 
The  irreversibility  condition  T(t)  C  T(f  +  At)  is  enforced  in  the  strong-form  equations  by  introducing  a 
strain-history  field,  H,  which  satisfies  the  Kuhn-Tucker  conditions  for  loading  and  unloading 

ipt-n<  o,  n>  o,  n(ipt-n)=  o  (in 


(see  Miehe,  Hofacker,  and  Welschinger  (2010a)  for  motivation  on  the  introduction  of  TL).  After  substituting 
'H  for  ipt  in  (9)2  we  get  the  modified  strong  form  equations  of  motion 


(SW 


/4e(l  -k)H 

V  0c 


day 
Ox  ; 


=  PUi 


2  d2C 

c~4edtf~ 


on  fix]0,T[ 
on  fix]0,T[. 


The  equations  of  motion  are  subject  to  the  boundary  conditions 


(12) 


(S:  BC)  ^ 


^ij^j  hi 
dc 


dxj 


rii  =  0 


on  dtig.x]0,T[ 
on  dfthi  x]0,T[ 

on  9(7x]0,T[ 


(13) 


with  gi(x ,  t)  and  hi(x ,  t)  being  prescribed  on  d£lgi  and  dflhi9  respectively,  for  all  t  e]0,  T[,  and  with  n(x) 
being  the  outward-pointing  normal  vector  of  the  boundary. 

In  addition,  the  equations  of  motion  (12)  are  supplemented  with  initial  conditions 


u(x,  0)  =  Uo(x) 

x  £  Cl 

ii{x ,  0)  =  vq(x ) 

x  e  Cl 

(14) 

c(x,  0)  =  co(a) 

X  £  Cl 

for  both  the  displacement  field  and  the  phase-field.  The  initial  phase-field,  co(x),  can  be  used  to  model 
pre-existing  cracks  or  geometrical  features  by  setting  it  locally  equal  to  zero  (see  Appendix  A  for  details). 


2.3  Analytical  solutions  of  the  one-dimensional  quasi-static  problem 

To  illustrate  various  properties  of  the  phase-field  formulation  for  brittle  fracture,  we  study  the  analytical 
solution  to  the  boundary  value  problem  introduced  in  the  previous  section.  We  restrict  ourselves  to  the  one¬ 
dimensional  domain  Cl  =  R  (d  =  1)  and  ignore  all  temporal  derivatives.  In  addition  we  assume  the  parameter 
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k  to  be  equal  to  zero,  and  the  strain  field  to  be  non-negative  (i.e.,  ipe  =  0  such  that  =  \ c2Ee 2).  Under 
these  assumptions,  the  strong  form  equilibrium  equations  (12)  reduce  to 


< 


c  —  4t 


d2^ 


dx2 


=  1 


on  Mx]0,T[ 
on  Rx]0,T[ 


with  a  =  c2Es.  By  virtue  of  (15)i,  the  stress  cr  is  constant  over  the  domain. 


(15) 


2.3.1  Homogeneous  solution 

We  first  study  the  homogeneous  solution  by  ignoring  all  spatial  derivatives  of  c.  From  the  phase-field  equation 
we  then  obtain 


Chom  —  \ 


(  2eE_  _2 
\Sc  • 


-1 


"hom 


+  1 


-1 


=  T-L  (loading) 

<  T-L  (unloading) 


(16) 


where  Chom  and  £hom  are  the  homogeneous  phase-field  and  strain,  respectively.  Substitution  of  this  result 
into  the  constitutive  equation  yields  the  homogeneous  stress  as  a  function  of  the  homogeneous  strain 


&  hom  —  \ 


\  -2 

J  Ee hom  (loading) 

[  (if  ^  +  l)  Ee\mm  ip+  <  H  (unloading) 


(  2eE  c2 

\~d7£h 


+  1 

-2 


(17) 


Figure  2(a)  shows  a  characteristic  plot  of  the  homogeneous  stress  versus  the  homogeneous  strain.  Figure  2(b) 
shows  the  corresponding  evolution  of  the  homogeneous  phase-field  parameter.  Note  that  the  plotted  results 
have  been  non-dimensionalized  as  outlined  in  Appendix  B.  It  is  observed  that  as  the  strain  is  increased, 
initially  the  stress  also  increases.  This  increase  in  stress  is  accompanied  by  a  gradual  decrease  of  the  phase- 
field  parameter.  At  some  point,  a  critical  stress  level,  ac,  is  reached  after  which  both  the  stress  and  the 
phase-field  decrease  in  value  upon  an  increase  in  strain.  In  unloading  the  phase-field  remains  constant,  which 
results  in  secant  unloading  behavior  in  the  stress-strain  curve. 

As  the  critical  stress  indicates  the  state  at  which  the  material  starts  to  soften,  we  refer  to  crc  as  the  crack 
nucleation  stress.  The  critical  value  for  the  stress,  and  corresponding  value  for  the  strain,  are  found  as 


9  [Wc 
16V“67’  £c 


(18) 


From  these  expressions  it  is  observed  that  the  crack  nucleation  stress  will  increase  as  e  decreases.  In  the 
limit  as  e  goes  to  zero,  i.e.,  when  the  phase-field  formulation  coincides  with  the  discrete  fracture  formulation, 
the  crack  nucleation  stress  becomes  infinite.  This  observation  is  consistent  with  the  properties  of  Griffith’s 
theory,  which  only  allows  for  crack  nucleation  at  stress  singularities.  It  is  interesting  to  note  that  the  critical 
value  for  the  phase-field  is  independent  of  the  model  and  material  parameters 


Cc  — 


3 

4 


(19) 


This  implies  that,  no  matter  what  parameters  are  used,  a  25  percent  reduction  in  the  phase-field  is  established 
prior  to  crack  nucleation  at  places  where  a  crack  emerges.  At  places  where  no  crack  is  formed,  the  value 
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(a)  (b) 

Figure  2:  One-dimensional  characteristic  stress-strain  (a)  and  phase- field- strain  (b)  curves  for  the  homoge¬ 
neous  solution.  Note  that  the  value  of  e  influences  the  maximum  tensile  stress,  see  (18). 


of  the  phase-field  scales  with  the  ratio  of  the  maximum  stress  in  that  point  and  the  crack  nucleation  stress. 
Consequently,  as  the  model  parameter  e  approaches  zero,  the  phase-field  will  approach  a  value  of  one  outside 
of  the  fracture  zone. 

By  virtue  of  the  preceding  analysis,  we  view  e  as  a  material  parameter  since  it  influences  the  critical  stress 
at  which  crack  nucleation  occures. 


2.3.2  Non-homogenous  solution 

Additional  interesting  features  of  the  model  can  be  observed  from  the  non-homogeneous  solution  of  the  one¬ 
dimensional  static  problem  (15).  If  we  ignore  the  irreversibility  in  the  model,  i.e.  we  take  H  =  \Ee2,  we  can 
combine  the  constitutive  and  phase-field  equations  to  yield  the  non-linear  ordinary  differential  equation 


/  2ecr2 
\c  4Egc 


2  d2c 

C  “  4C  “  ' 


(20) 


A  solution  for  the  phase-field,  c(x)9  with  a  crack  at  x  =  0  is  found  by  supplementing  this  differential  equation 
with  the  boundary  condition  lim  c(x)  =  Chom(^),  which  implies  that  the  phase-field  gradient  vanishes  far 

£—>•±00 

away  from  the  crack.  It  is  important  to  note  that  the  homogeneous  solution  to  the  phase-field  is  dependent  on 
the  stress,  cr,  as  elaborated  in  the  previous  section.  In  addition  to  the  far-field  boundary  condition,  we  require 
the  solution  to  be  symmetric  and  differentiable  at  every  point  except  for  x  =  0. 

The  first  step  in  finding  the  non-homogeneous  solution  to  the  phase-field  problem  is  to  multiply  equation 
(20)  with  ^  and  make  use  of  the  fact  that,  by  (15),  cr  is  constant,  to  obtain 


d 

ecr2  c 2  2  | 

fdc\2 

dx 

£EQC  +  2  26  * 

VdaJ 

(21) 


Since  we  require  the  solution  to  be  symmetric  around  x  =  0,  we  integrate  this  expression  from  x  to  infinity 
for  positive  values  of  x9  and  from  minus  infinity  to  x  when  x  is  negative.  Since  we  have  specified  the  crack 
to  be  centered  at  x  =  0,  we  require  the  phase-field  to  have  a  minimum  at  x  =  0.  From  these  requirements, 
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we  obtain 


dc 

dx 


=  sgn(x) 


ecr* 


c2Egc 


—  c  —  a 


where  the  coefficient  a  follows  from  the  far-field  boundary  condition  as 


6(7 


a  =  —  - 


uhom 


egc 


+ 


uhom 


^hom  • 


(22) 


(23) 


By  definition,  substitution  of  the  homogeneous  solution  c^om(cr)  into  equation  (22)  yields  a  zero  phase- 
field  gradient.  An  analytical  non-homogeous  solution  to  equation  (20)  can  be  found  for  the  case  of  a  fully 
developed  crack,  i.e.  for  a  =  0  (no  stress  carrying  capability)  and  Chom(0)  =  1.  In  this  case,  equation  (22) 
reduces  to 

dc  .  N 1  —  c 

S  =  8gn(x>^T'  <24) 

and,  assuming  c(0)  =  0,  we  find 

c=l-exp^-^^  (25) 

as  the  solution  that  satisfies  the  specified  boundary  conditions.  Note  that  this  is  the  exact  function  used  by 
Miehe  et  al.  (2010a)  to  derive  their  phase  field  formulation. 

The  unique  non-homogeneous  solution  can  also  be  constructed  for  a  >  0.  In  this  case,  a  second  admiss- 
able  phase-field  value,  smaller  than  the  homogeneous  solution,  Chom,  for  which  the  gradient  is  equal  to  zero 
can  be  found.  This  phase-field  value  corresponds  to  the  value  of  the  phase-field  at  the  center  of  the  crack, 
and  is  denoted  by  ccrack.  Since  the  phase-field  increases  monotonically  from  ccrack  at  the  center  of  the  crack 
(at  x  =  0)  to  Chom  far  away  from  the  crack  (x  =  ±oo),  we  find  the  coordinates  0  <  ±x  <  oo  at  which  the 
phase-field  is  equal  to  ccrack  <  c(x)  <  Chom  by  evaluation  of 


c(x) 


x  =  =b 


c2 


1  (  €<J 

2?  \WWc  +  "2  ~c~a 


dc. 


(26) 


Ccrack 


In  non-dimensional  form  with  the  non-dimensional  constant  Ch  =  1,  (26)  becomes 

c(L0x *) 


*  I 

X  =  ± 


/ 


l 


2(e*)2 


e*(a*)2  .  c 


+  --c-a 


dc. 


(27) 


Evaluating  this  integral  numerically  for  a  =  /3  a  c  with  various  values  of/?  E]0, 1[  we  get  the  solutions  shown 
in  Figure  3  (note  that  the  presented  results  have  again  been  non-dimensionalized  and  merely  depend  on  the 
parameter  /?).  It  is  interesting  to  note  that,  except  for  the  limiting  case  of  a  fully  developed  crack,  the  solution 
to  the  phase-field  is  smooth  at  the  center  of  the  crack. 


3  Numerical  formulation 

The  numerical  solution  of  (12)  requires  a  spatial  and  temporal  discretization.  In  this  section  we  formulate 
the  spatial  discretization  by  means  of  the  Galerkin  method  and  we  introduce  two  temporal  discretization 
schemes:  a  monolithic  implicit  scheme  and  a  staggered  scheme  in  which  the  momentum  equation  and  phase- 
field  equation  are  solve  separately,  and  in  which  the  momentum  equation  may  be  solved  either  by  an  explicit 
or  implicit  scheme. 
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Figure  3:  One-dimensional  solution  of  the  phase  field  formulation  for  various  values  of  the  stress  ratio  (3  = 
—  =  ^r.  Note  that,  except  for  the  limiting  case  (3  =  0,  the  phase-field  is  smooth  at  the  center  of  the  crack 

(x*  =  0). 


3.1  Continuous  problem  in  the  weak  form 

For  the  weak  form  of  the  problem  we  define  the  trial  solution  spaces  St  for  the  displacements  and  St  for  the 
phase-field  as 

St  =  {u(t)  e  (Ff1(n))d  |  Ui(t)  =  gt  on  dQgi}  (28) 

St  =  {c(t)  €  H\Q)}.  (29) 

Similarly,  the  weighting  function  spaces  are  defined  as 

V  =  {w  G  |  Wi  =  0  on  dCl9i}  (30) 

V  =  {qeH1(Q)}.  (31) 

Multiplying  the  equations  in  (12)  by  the  appropriate  weighting  functions  and  applying  integration  by  parts 
leads  to  the  weak  formulation: 


Given  g ,  h ,  no,  no,  and  Co  find  u(t)  G  St  and  c(t)  G  St,  t  G  [0,  T],  such  that  for  all 
w  G  V  and  for  all  q  G  V, 


(WW 


f  Ue{l-k)n 

Vv 


{pit,  w)n  +  (<T,  Vw)fi  =  (h,  w)dQh 


+  (4e2Vc,  =  (l,g)n 


(pu(0),w)n  =  (pu0,w)n 
(pu(0),w)n  =  (pu0,w)n 
(c(0  ),q)n  =  (co,q)a 


(32) 


where  (•,  «)n  is  the  £2  inner  product  on  Cl. 


3.2  The  semidiscrete  Galerkin  form 

Following  the  Galerkin  method,  we  let  S ^  C  St,  Vh  C  V,  Sj1  C  5t,  and  Vh  C  V  be  the  usual  finite¬ 
dimensional  approximations  to  the  function  spaces  of  the  weak  form  (see  Hughes  (2000)  for  details).  The 
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semidiscrete  Galerkin  form  of  the  problem  is  then  given  as 

Given  g ,  h ,  no,  Uo,  and  Co  find  uh(t )  G  and  c^(t)  G  <S/\  t  G  [0,  T],  such  that  for 
all  wh  e  Vh  and  for  all  qh  eVh, 

(puh,wh^  +  (cr,\7wh)Q  =  (h,wh)dQh 
((4e(1-fc)^  +  ^  ^ ^ ^  +  (4e2 Vcfc,  v^)n  =  (1, (33) 

(puh(0),wh)Q  =  ( pu0,  wh)n 
(puh(0),wh^  =  ( pii0,wh)n 
(ch(0),qh)Q  =  (c0,  qh)n 

The  explicit  representations  of  uh ,  c^,  and  qh  in  terms  of  the  basis  functions  and  nodal  variables  are 


nb 


Ui  =  ^NA(x)diA 

A 

nb 

(34) 

wi  =  ^ ^NA(x)ciA 

A 

nb 

(35) 

Ch  =  '^jNA(x)(pA 

A 

nb 

(36) 

qh  =  ^2 Na(x)xa 

(37) 

A 


where  n 5  is  the  dimension  of  the  discrete  space,  the  Na  s  are  the  global  basis  functions,  i  is  the  spatial 
degree-of-freedom  number,  and  diA,  cia ,  4>a,  and  xa  are  control  variable  degrees-of-freedom.  Note  that  we 
have  assumed  that  both  the  finite  dimensional  trial  solution  and  weighting  function  spaces  are  defined  by  the 
same  set  of  basis  functions. 

3.2.1  Isogeometric  spatial  discretization 

In  contrast  to  earlier  work  on  phase-field  models,  we  have  chosen  to  use  isogeometric  spatial  discretizations  as 
introduced  by  Hughes  et  al.  (2005),  which  are  based  on  NURBS  and  T-splines.  T-splines  are  a  generalization 
of  NURBS  that  allows  greater  flexibility  in  geometric  design  including  local  refinement.  A  class  of  analysis- 
suitable  T-splines  was  identified  in  Li  et  al.  (2010)  and  Scott  et  al.  (201 1).  This  class  of  T-splines  preserves  the 
important  mathematical  properties  of  NURBS  while  providing  an  efficient  and  highly  localized  refinement 
capability.  Analysis-suitable  T-splines  possess  the  following  properties: 

•  Linear  independence. 

•  The  basis  constitutes  a  partition  of  unity. 

•  Each  basis  function  is  non-negative. 

•  An  affine  transformation  of  an  analysis- suitable  T-spline  is  obtained  by  applying  the  transformation  to 
the  control  points.  We  refer  to  this  as  affine  covariance.  This  implies  that  all  patch  tests  (see  Hughes 
(2000))  are  satisfied  a  priori. 
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•  They  obey  the  convex  hull  property,  see  Piegl  and  Tiller  (1997). 

•  Local  refinement  is  possible. 

For  additional  details  on  basic  T-spline  analysis  technology  see  Scott  et  al.  (2010). 

The  spline-based  analysis  strategy  of  isogeometric  analysis  has  shown  some  advantages  when  compared 
to  standard  C°  finite  elements.  First,  isogeometric  analysis  allows  for  the  efficient  and  exact  geometric  rep¬ 
resentation  of  many  objects  of  engineering  interest.  Also,  when  using  T- spline-based  isogeometric  analysis, 
efficient  and  automatic  local  mesh  refinement  strategies  exist  that  preserve  the  exact  geometry.  Secondly, 
the  isogeometric  basis  is  in  general  smooth.  Although  the  phase-field  model  permits  the  use  of  traditional 
C°  finite  elements,  the  use  of  a  smooth  base  is  anticipated  to  have  favorable  effects.  One  is  that  stresses 
are  represented  more  accurately  than  with  traditional  C°  finite  elements,  which  has  been  observed  to  yield 
efficient  spatial  discretizations  for  discrete  (see  Verhoosel  et  al.  (2010))  and  smeared  (see  Verhoosel  et  al. 
(2011))  fracture  models.  Another  study  by  Benson  et  al.  (2011)  has  shown  that  the  accuracy  of  the  spatial 
discretization  of  NURBS  allows  for  more  efficient  time  integration  for  explicit  dynamics  in  the  context  of 
large  deformation  shells. 


3.3  Time  discretization  and  numerical  implementation 

3.3.1  Monolithic  generalized-^  time  discretization 


The  monolithic  time  integration  scheme  is  based  on  the  generalized-^  method  introduced  by  Chung  and 
Hulbert  (1993).  We  define  the  residual  vectors  as 

R“  =  {RUA,i},  (38) 

R\i  =  (h,  NAei)dUh  -  ( pu\  NAe., i)fi  -  (ajk,  Bf  ^  (39) 

and 


Rc  =  {Rca}, 

RcA  =  (l,NA)n-^ 


4e(l  -  k)U 


+  l)c'I,iVA)  - 
n 


where  e,  is  the  ith  Euclidean  basis  vector  and 


bT  - 


=  1  (dNA 

2  V  dxj 


r  ,9NA 

Oik  +  - Oi 

dxk  ■ 


/  3Na\ 

V  6  dxS  dxi 


(40) 

(41) 


(42) 


so  that  £jk  =  BlikdA,i.  For  time  step  n,  let  dn  and  (j)n  be  the  vectors  of  control  variable  degrees-of- 
freedom  of  the  displacements  and  phase-field,  respectively  (see  (34)  and  (36)).  We  then  define  vn  =  dn,  and 
an  =  dn.  The  monolithic  generalized- a  time  integration  scheme  is  then  stated  as  follows:  given  (dn,  vn, 


find  (dn+i,  vn+i,  an- 1_  1 ,  dn-|_a^,  Vn+a/,  o-n+am,  0n+i)  such  that 

FI  (dn_|_a^ ,  vn-\-aj ,  <2n+am ,  0n_|_^)  =  0,  (43) 

F lc(dn+a/ ,  </>n+1)  =  0,  (44) 

dn-\-otf  —  dn  T"  (Xfi^dn- |_i  dn),  (45) 

Vn+cxf  =Vn  +  af(vn+1  -  Vn),  (46) 

^n+am  =  Un  +  CVm(<2n_|_i  O'n)')  (47) 

Vn+ 1  =  Vn  +  Af((l  -  7 )an  +  7an+i),  (48) 

(At)2 

dn+i  =  dn  +  A tvn  4 - — — ((1  —  2 f3)an  +  2/4an+i),  (49) 
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where  At  =  tn+ 1  —  tn  is  the  time  step  and  the  parameters  af,am,l3,  and  7  define  the  method.  These 
parameters  will  be  discussed  below. 

At  each  time  step,  the  solution  is  obtained  using  a  Newton-Raphson  method  to  solve  the  nonlinear  equa¬ 
tions  above.  Letting  i  be  the  Newton  iteration,  the  residual  vector  and  consistent  tangent  matrix  for  the 
linearized  system  are  defined  by 


where 


ony 

ddn-\- 1 
<9Ri 

£^n+l 


Aoj  T~ 

A  CL  ~\~ 


dny 

d(t>n+ 1 

d'fin+l 


A</>  =  -R", 
A  0  =  -Rf, 


R“  =  R“(d; 


n+a  f >  ^n+a  f  i 


>  ^n+l)? 


and 


(50) 

(51) 


(52) 


R;  =R  c(dn+a/,0n+1).  (53) 

For  each  iteration,  the  linearized  system  defined  by  (50)  and  (51)  is  solved  and  iteration  continues  until 
convergence  of  the  residual  vectors  occurs.  For  the  examples  discussed  below,  we  have  defined  convergence 
as 


max 


iiRrii  iiRiii 


<  tol 


(54) 


where  ||  •  ||  denotes  the  Euclidean  norm.  For  the  simulations  considered  in  this  work,  tol  =  10  4  has  been 
observed  to  be  an  appropriate  choice. 


Parameter  selection  In  Chung  and  Hulbert  (1993)  it  was  shown  that  af  and  am  can  be  parametrized 
by  the  spectral  radius,  p oo,  of  the  amplification  matrix  at  Af  =  00  such  that  second-order  accuracy  and 
unconditional  stability  are  achieved  for  a  second-order  linear  problem  if 


°if 


C^rri 


1 

Poo  +  1  ’ 
2  Poo 
Poo  +  1  ’ 


P=  |(l+£»m-a/)2, 

1 

7  —  2  • 


(55) 

(56) 

(57) 

(58) 


3.3.2  Staggered  time  discretization 

For  the  staggered  time  integration  scheme,  the  momentum  and  phase-field  equations  are  solved  indepen¬ 
dently.  At  a  given  time  step,  the  momentum  equation  is  solved  first  to  get  the  displacements.  Using  the 
updated  displacements,  the  phase-field  equation  is  solved.  In  addition  to  reducing  the  problem  to  solving  two 
linear  systems,  this  scheme  also  allows  greater  flexibility  in  how  the  momentum  equation  is  solved,  i.e.,  we 
can  use  either  implicit  or  explicit  schemes.  This  scheme  can  also  be  generalize  to  a  predictor/multicorrector 
format  where  additional  Newton-Raphson  iterations  can  be  performed  within  a  time  step.  Below  we  present 
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a  general  predictor/multicorrector  algorithm,  but  for  the  results  presented  later  we  use  only  one  pass  of  the 
corrector  stage. 

Defining  the  residual  vectors  for  the  momentum  and  phase-field  equations  by  (39)  and  (41),  and  again  let¬ 
ting  d  and  <j>  be  arrays  of  the  control  variable  coefficients  in  (34)  and  (36),  the  staggered  predictor/multicorrector 
time  integration  scheme  is  stated  as  follows:  given  (dn,  vn,  an,  0n),  solve 

Predictor  stage 


Multicorrector  stage 


i  —  0  (iteration  counter) 

Vn+1  =vn  +  At(l  -  7 )an 
(At)2 


dn- (-1  dn 


A  tvr 


°n+l  =  0 

(i) 

vn+l  =  ^n+l 


°^n+ 1  —  dn+ 1 


/W 

ln+ 

=  <t>n 


(1  -  2 P)an 


(an+l  an) 


(i) 

an+am  ~  an  ~ f 
Vn+a/  =Vn  +  af(vn+l  ~  Vn ) 

dn+ct  f  dn  +  L 


dn) 


M*Aa  =  R  u(d[ 


a 


(i+l) 

n+1 


=  a 


« 

'n+1 


n(0 

n+a  f  ^n+a /  )  u,n+a: 

Aa 


(i+1)  ~  ,  a  -L  (i+1) 

<+i  =  vn+i  +  At7a„+i 

dn+1  =  dn+ 1  +  (A+’.fa+j0 

KccA</>  =  Fc 

W 


,<t>[ 


(*)  \ 
n+ll 


(59) 

(60) 

(61) 

(62) 

(63) 

(64) 

(65) 


(66) 

(67) 

(68) 

(69) 

(70) 

(71) 

(72) 

(73) 

(74) 


The  phase-field  arrays  in  (73)  are  defined  as 


Kcc 

Kab 


[Kab], 


//4e(l  -k)H 
U  Qc 


Fc  =  {Fa}, 
Fa  =  (1,Na), 


Nb,Na 


ONa 

dxi 


Q 


(75) 

(76) 

(77) 

(78) 


and  At  =  tn+ i  —  tn  is  the  time  step  and  the  parameters  am,af,(3,  and  7,  which  define  the  method,  are 
selected  as  described  below. 

If  the  linearized  momentum  equation  (69)  is  being  solved  implicitly  then 

F) r>  u 

M*  =  —  — — —  =  amM  +  a//?(Ai)2K,  (79) 

dan+ 1 
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where  M  is  the  consistent  mass  matrix  and 

K  =  [KYb,hI  (80) 

KAB,iJ=(^rLBBmn^k)  (81) 

is  the  consistent  damage-elastic  tangent  stiffness  matrix.  If  we  let  M*  =  amM.  where  M  is  the  lumped 
mass  matrix  then  the  linearized  momentum  equation  is  solved  explicitly.  When  computing  with  NURBS  and 
T-splines,  we  compute  M  using  the  row-sum  technique  as  described  in  Hughes  (2000).  Due  to  the  fact  that 
NURBS  and  T- spline  basis  functions  are  point- wise  positive,  the  row-sum  lumped  mass  matrix  is  guaranteed 
to  be  positive.  Furthermore,  it  is  also  mass  conservative. 


Parameter  selection  For  the  staggered  solution  strategy,  the  choice  of  am,  af ,  and  M*  provides  several 
options  for  the  type  of  algorithm  that  is  used  to  solve  the  linear  momentum  problem.  For  the  fully  implicit  case 
(M*  =  amM  +  af/3(At)2K.)  we  use  the  generalized-^  method  described  above.  For  the  fully  explicit  case 
(M*  =  amM)  we  use  either  the  HHT-a  of  Hilber,  Hughes,  and  Tayler  (1977)  or  the  explicit  generalized-^ 
method  of  Hulbert  and  Chung  (1996).  The  HHT-a  method,  parameterized  by  a ,  provides  a  second-order 
accurate  family  of  algorithms  for  linear  second-order  equations  if  a  G  0]  and 


m  — 

i, 

(82) 

lf  = 

1  +  OL-, 

(83) 

(1  —  a)2 

P  = 

4 

(84) 

1  -  2a 

7  = 

2 

(85) 

(see  Miranda,  Ferencz,  and  Hughes  (1989)).  The  explicit  generalized-a  method,  as  shown  by  Hulbert  and 
Chung  (1996),  is  a  one-parameter  family  of  explicit  algorithms  that  provides  optimal  numerical  dissipation 
and  is  second-order  accurate  for  linear  problems  if  a f  =  0  and 


(%rn 

(3 


7 


2  -  pb 

1  +  Pb  ’ 

5  -  3  pb 

(l  +  p6)2(2-p6)’ 

1 

am  +  2> 


(86) 

(87) 

(88) 


where  pi,  is  the  spectral  radius  value  at  the  bifurcation  limit  of  the  principal  roots  of  the  characteristic  equation. 


4  Numerical  results 

In  this  section  we  investigate  the  numerical  performance  of  the  phase-field  fracture  model.  All  geometries 
have  been  discretized  spatially  using  either  NURBS  or  T- spline  basis  functions  (which  we  refer  to  as  the 
global  smooth  basis).  The  numerical  computation  for  all  the  models  was  performed  using  the  Bezier  ex¬ 
traction  methods  described  by  Borden  et  al.  (2010)  and  Scott  et  al.  (2010).  Bezier  extraction  constructs  the 
minimal  set  of  Bezier  elements  defining  a  NURBS  or  T- spline.  A  Bezier  element  is  a  region  of  the  physical 
domain  in  which  the  basis  functions  are  C°° -continuous  and  over  which  integration  is  performed.  In  addition 
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to  defining  elements,  Bezier  extraction  also  builds  an  extraction  operator  for  each  Bezier  element  that  maps 
a  Bernstein  polynomial  basis  defined  on  the  Bezier  element  to  the  global  smooth  basis.  The  transpose  of  the 
extraction  operator  maps  the  control  points  of  the  global  smooth  basis  to  the  Bezier  control  points.  The  idea 
is  illustrated  in  Figure  4  for  a  cubic  B-spline  curve.  The  B-spline  curve  with  its  control  points,  P  =  {P /}/=1, 
and  basis  functions,  N  =  {N/}j=1,  is  shown  on  the  left  with  the  elements,  Ue,  defined  on  the  parametric  do¬ 
main.  On  the  right  we  illustrate  the  action  of  the  extraction  operator,  C2  G  M4x4,  for  U2.  The  transpose  of  the 
extraction  operator,  C^,  defines  the  Bezier  control  points,  Q2  =  {Q2j}j=1,  for  the  Bezier  element  in  terms 
of  the  global  control  points,  as  shown  in  Figure  4.  The  Bernstein  polynomial  basis  functions,  B  =  {B i}j=1, 
are  shown  below  the  Bezier  element.  The  extraction  operator  is  used  to  extract  the  smooth  B-spline  basis 
functions  N2  =  {N/}f=2  from  the  Bernstein  basis,  i.e.,  the  B-spline  basis  functions  can  be  written  as  linear 
combinations  of  the  Bernstein  basis  functions.  Bezier  extraction  provides  a  completely  local  representation 
of  the  global  smooth  basis.  It  provides  an  element  data  structure  that  can  be  integrated  into  existing  finite 
element  frameworks  in  a  straightforward  manner. 

To  integrate  arrays  over  the  Bezier  elements  we  use  Gaussian  quadrature  with  a  p  +  1  rule  in  each  para¬ 
metric  direction,  where  p  is  the  polynomial  degree  of  the  basis  functions.  Thus,  in  two-dimensions  we  use 
a  3 -by- 3  quadrature  rule  for  quadratic  basis  functions  and  a  4-by-4  quadrature  rule  for  cubic  basis  functions. 
Hughes,  Reali,  and  Sangalli  (2010)  have  shown  that  smooth  bases  allow  for  more  efficient  quadrature  rules. 
See  also  the  appendix  in  Hughes,  Reali,  and  Sangalli  (2008)  which  demonstrates  the  stability  and  accuracy  of 
reduced  quadrature  rules  for  NURBS.  This  is  beyond  the  scope  of  the  work  presented  here  and  we  acknowl¬ 
edge  that  the  quadrature  rules  we  use  may  represent  overkill. 

For  the  examples  below,  the  reported  mesh  sizes,  h ,  are  computed  on  the  Bezier  elements  as  h  =  tfa 
where  a  is  the  area  of  an  element  in  two  dimensions  and  the  volume  of  an  element  in  three  dimensions  and  d 
is  the  number  of  spatial  dimensions.  In  most  cases,  the  mesh  is  such  that  h  =  e/2  in  the  area  where  a  crack 
has  formed.  Experience  has  shown  that  this  relationship  between  h  and  e  provides  sufficient  accuracy  without 
over  resolving  the  crack. 

4.1  2D  Quasi-static  shear  load 

In  this  section  we  consider  a  quasi-static  benchmark  test  from  Miehe,  Hofacker,  and  Welschinger  (2010a) 
in  order  to  compare  the  results  obtained  from  standard  C°  finite  elements  and  isogeometric  finite  elements. 
The  geometry  and  boundary  conditions  of  the  model  are  shown  in  Figure  5.  We  use  C2 -continuous  cubic 
T-splines  for  the  spatial  discretization  of  the  model  and  the  staggered  quasi-static  solution  strategy  described 
by  Miehe,  Hofacker,  and  Welschinger  (2010a)  to  obtain  the  solution  at  each  load  increment.  As  can  be  seen 
from  Figure  5b,  T-splines  can  be  locally  refined  in  the  area  where  the  crack  forms.  This  is  in  contrast  to  a 
NURBS,  where  refinement  propagates  globally  to  produce  a  much  denser  mesh. 

The  material  parameters  are  E  =  210  GPa,  v  =  0.3,  and  Qc  =  2700  J/m2  and  plane  strain  is  assumed. 
The  length  scale  is  chosen  to  be  e  —  7.5  x  10-6m  and  we  do  not  include  a  viscous  damping  term  on  the 
phase-field.  The  T-spline  was  refined  a  priori  based  on  the  expected  solution.  The  initial  crack  is  modeled  as 
a  discrete  discontinuity  in  the  geometry  and  the  T-spline  contains  C°  lines  that  radiate  out  from  the  crack  tip 
so  that  the  mesh  is  divided  into  four  equal  square  subdomains  of  C2  continuity.  These  C°  lines  were  included 
to  facilitate  modeling  the  sharp  crack  tip  as  a  C°  geometric  feature.  An  alternative  would  have  been  to  use 
a  globally  C2 -continuous  T-spline  on  the  entire  domain  and  to  introduce  the  crack  through  an  induced  crack 
in  the  phase-field  (see  Appendix  A).  The  locally  refined  T-spline  contained  5,587  cubic  basis  functions.  The 
calculations  of  Miehe,  Hofacker,  and  Welschinger  (2010a)  utilized  30,000  linear  triangles. 

Figure  6  shows  the  progression  of  the  crack  at  several  load  levels  and  the  load-displacement  curve  is 
shown  in  Figure  7  with  a  comparison  to  the  results  reported  by  Miehe,  Hofacker,  and  Welschinger  (2010a). 
As  can  be  seen  from  the  load-displacement  curve,  the  results  obtained  from  the  T-spline  mesh  are  in  good 
agreement  with  those  obtained  from  standard  C°  finite  elements,  but  with  far  fewer  degrees -of-freedom. 
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Physical  Space 


P5 


Global  view 


Element  view 


Figure  4:  Bezier  extraction  for  a  cubic  B- spline  curve.  The  B- spline  curve  and  basis  functions  are  shown  on 
the  left.  The  action  of  the  extraction  operator,  C2,  for  element  ^2  is  illustrated  on  the  right.  The  transpose  of 
the  extraction  operator  defines  the  control  points,  Q2  =  {Q2,/}/=i>  °f  the  Bezier  element  from  the  control 
points,  P2  =  {P/}/=2  °f  the  B-spline  curve.  The  B-spline  basis  functions,  N2  =  {N/}f=2,  can  be  computed 
over  the  element  by  applying  the  extraction  operator  to  Bernstein  polynomial  basis  functions,  B  =  {B 
defined  on  the  Bezier  element.  Note  that  the  Bernstein  basis  is  the  same  for  each  element.  Formation  of 
element  arrays  can  thus  be  standardized  in  a  shape  function  routine. 


We  attribute  this  to  the  smooth  description  of  the  stress  fields,  which  was  also  observed  to  be  beneficial  in 
cohesive  zone  modeling  by  Verhoosel  et  al.  (2010)  and  gradient  damage  modeling  by  Verhoosel  et  al.  (201 1). 
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Figure  5 :  Input  model  for  the  quasi-static  shear  benchmark  test:  (a)  geometry  and  boundary  conditions  and 
(b)  the  Bezier  element  representation  of  the  T- spline.  The  T- spline  contains  5587  cubic  basis  functions  and  the 
effective  element  size  is  between  hmin  =  3.906  x  10-3  mm  and  /imacc  =  6.25  x  10-2  mm.  The  refinement 
was  performed  a  priori  based  on  knowledge  of  the  expected  crack  path. 


u  =  9.0  x  10  3  mm  u  =  11.0  x  10  3  mm  u  =  13.4  x  10  3  mm 

Figure  6:  Crack  progression  for  the  quasi-static  shear  test  at  several  load  increments. 
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Displacement  u  [mm] 


Figure  7:  Force-displacement  curve  for  the  quasi-static  shear  test.  The  results  obtained  using  the  T-spline 
discretization  are  compared  to  those  reported  by  Miehe,  Hofacker,  and  Welschinger  (2010a)  using  standard 
C°  finite  elements. 
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4.2  Dynamic  crack  branching 

In  this  example,  we  model  a  pre-notched  rectangular  plate  loaded  dynamically  in  tension.  The  geometry 
and  boundary  conditions  of  the  problem  are  shown  in  Figure  8.  A  traction  load  is  applied  to  the  top  and 
bottom  surface  at  the  initial  time  step  and  held  constant  throughout  the  simulation.  All  other  surfaces  have 
a  zero  traction  condition  applied.  This  load  condition  is  such  that  crack  branching  will  occur  (see  Song, 
Wang,  and  Belytschko  (2008)  for  a  report  of  results  for  this  problem  using  several  other  methods  of  dynamic 
fracture  analysis).  The  initial  crack  is  induced  by  an  initial  strain-history  field  (see  Appendix  A)  allowing  the 
geometry  to  be  modeled  as  a  continuous  quadratic  C1  -continuous  NURBS  patch. 
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Figure  8:  The  geometry  and  boundary  conditions  for  the  crack  branching  example.  In  the  simulation  the 
initial  crack  is  modeled  by  introducing  an  initial  strain  history  field  that  induces  a  phase  field  at  the  initial 
crack  location,  and  the  geometry  and  displacement  field  are  C 1  continuous  throughout  (see  Appendix  A). 

The  model  parameters  are  p  =  2450  kg/m3,  E  =  32  GPa,  v  =  0.2,  Qc  =  3  J/m2,  and  plane  strain 
is  assumed.  The  corresponding  dilatational,  shear,  and  Rayleigh  wave  speeds  are  Vd  =  3810  m/s,  vs  = 
2333  m/s,  vr  =  2125  m/s.  The  length  scale  was  chosen  to  be  e  =  2.5  x  10-4  m.  The  monolithic  generalized- 
a  time  integration  scheme  discussed  in  section  3.3.1  was  used  with  p ^  =  0.5. 

Table  1  lists  the  parameters  for  three  successively  finer  meshes  used  for  this  example.  A  uniform  mesh 
was  used  to  remove  any  effect  from  mesh  distribution  and  size  variation.  To  ensure  accurate  results,  the  time 
step  was  chosen  for  each  mesh  such  that  At  «  H/vr. 


h 

At 

Mesh  1 

2.5  x  10“4  m 

i  x  icr7  s 

Mesh  2 

1.25  x  10“4  m 

5  x  1CT8  s 

Mesh  3 

6.25  x  HT5  m 

2.5  x  1CT8  s 

Table  1 :  The  mesh  sizes  and  time  steps  used  for  the  dynamic  crack  branching  example.  The  mesh  is  a  uniform, 
quadratic  NURBS  in  each  case.  To  ensure  accurate  results,  the  time  step  size  is  chosen  to  be  roughly  h/vR. 

In  order  to  make  a  comparison  between  the  results  from  the  three  meshes  we  define  the  elastic  strain 
energy  as 

£e  -  /  {[(l  -  k)c2  +  k]ip+ (e)  +  tp~ (e)}  dx  (89) 

Jn 
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and  the  dissipated  energy  as 


£d  = 


( c  —  l)2  dc  dc 

4e  dxi  dxi 


dx. 


(90) 


These  values  provide  a  measure  of  the  global  response  of  the  model. 

The  final  phase-field  results  for  each  mesh  are  shown  in  Figure  9.  As  can  be  seen,  the  crack  path  is  similar 
for  all  meshes  although  mesh  1  (h  =  e)  fails  to  reach  the  boundary  in  the  same  amount  of  time  as  the  other 
two.  This  is  typical  of  the  behavior  seen  when  the  mesh  is  too  coarse  to  adequately  resolve  the  phase-field 
transition  near  the  crack.  Coarse  meshes  tend  to  limit  the  crack  tip  velocity.  It  is  difficult  to  see  a  difference 
in  the  crack  path  between  mesh  2  and  mesh  3. 


(a)  Mesh  1  with  h  =  2.5  x  10  4  m 


(b)  Mesh  2  with  h  =  1.25  x  10  4  m 


(c)  Mesh  3  with  h  =  6.25  x  10  5  m 

Figure  9:  Results  for  the  crack  branching  example  at  t  =  80  /is.  For  all  three  meshes  e  =  2.5  x  10-4  m. 

Figure  10(a)  shows  the  elastic  strain  energy  as  defined  by  (89)  and  Figure  10(b)  show  the  energy  dissipated 
by  the  formation  of  the  crack  as  defined  by  (90).  These  plots  show  that  h  =  e/ 2  provides  enough  resolution 
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to  capture  the  response  of  the  model.  For  h  =  e  the  mesh  is  too  coarse  and  the  material  behaves  as  if  it 
is  too  stiff  and  too  much  energy  is  dissipated.  This  is  a  result  of  the  mesh  being  too  coarse  to  capture  the 
gradients  of  the  phase-field  near  the  crack.  The  crack  tip  velocity  is  plotted  in  Figure  10(c).  In  all  cases,  the 
velocity  stays  well  below  60%  of  the  Rayleigh  wave  speed  as  has  been  commonly  observed  in  experiments 
(see  Ravi-Chandar  and  Knauss  (1984)). 
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Figure  10:  Plots  over  time  of  (a)  the  elastic  strain  energy  as  defined  by  (89),  (b)  the  dissipated  fracture  energy 
as  defined  by  (90),  and  (c)  crack  tip  velocity.  After  branching,  the  reported  crack  tip  velocity  is  that  of  the 
upper  branch. 

In  Figure  11  we  show  a  post-processed  plot  of  mesh  2  at  t  =  70  /is.  In  this  figure  we  have  scaled  the 
displacements  by  a  factor  of  50  and  removed  areas  of  the  model  from  the  plot  where  c  <  0.05  in  order  to 
show  a  representation  of  the  cracked  geometry. 


Remark:  For  the  monolithic  time  integration  scheme  the  size  of  At  is  limited  by  accuracy  and  convergence. 
For  the  dynamic  crack  branching  example  discussed  here,  choosing  a  time  step  such  that  At  <  2 Ii/vr 
produced  acceptable  accuracy.  For  At  =  SH/vr  the  Newton-Raphson  method  failed  to  converge  near  the 
time  of  the  initial  crack  propagation. 
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Figure  1 1:  A  post-processed  plot  of  mesh  2  at  t  =  70  /is.  The  displacements  have  been  scaled  by  a  factor  of 
50  and  areas  of  model  where  c  <  0.05  have  been  removed  from  the  plot.  Pressure  is  measured  in  Pascals. 


Remark:  Since  the  crack  is  not  tracked  algorithmically,  the  velocity  of  the  crack  tip  is  measured  as  a  post¬ 
processing  step.  For  the  values  reported  here,  the  location  of  the  crack  tip,  x ,  is  found  on  an  iso-curve  of  the 
phase-field  with  value  0.25.  The  velocity  is  then  computed  as  vn  =  (xn+i  —  xn)/At. 
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4.3  Dynamic  shear  loading 

In  this  example  we  model  crack  initiation  and  propagation  under  a  dynamic  shear  load.  The  model  is  based 
on  experimental  results  reported  by  Kalthoff  and  Winkler  (1987)  and  Kalthoff  (2000).  Previous  numerical 
results  of  this  problem  based  on  XFEM  have  been  reported  by  Belytschko  et  al.  (2003),  among  others,  and  a 
comparison  of  results  between  XFEM,  the  element  deletion  method,  and  the  interelement  crack  method  have 
been  reported  by  Song,  Wang,  and  Belytschko  (2008).  Numerical  results  of  a  similar  experiment  reported 
by  Ravi-Chandar  (1998)  have  been  reported  by  Remmers,  de  Borst,  and  Needleman  (2008)  where  cohesive 
segments  have  been  used  to  model  the  crack. 

The  input  geometry  and  loading  conditions  for  the  simulation  are  shown  in  Figure  12,  where  symmetry 
is  employed  to  reduce  the  computational  cost.  In  the  experiment,  the  load  was  applied  by  firing  a  projectile 
at  a  prenotched  specimen.  In  our  simulation  we  model  the  case  where  the  projectile  was  fired  with  a  velocity 
of  33  m/s  by  applying  the  kinematic  velocity 


v  = 


T0V o  t  <  ^ 


vo 


t>  to 


(91) 


where  vq  =  16.5  m/s  and  to  =  1ms-  A  no  traction  boundary  condition  is  applied  to  all  unspecified  sur¬ 
faces.  We  model  the  geometry  using  quadratic  NURBS  basis  functions.  The  initial  crack  is  modeled  by  a 
discontinuity  in  the  geometry  in  order  to  introduce  a  sharp  crack  tip  as  in  Section  4.1. 


100  mm 


Figure  12:  The  geometry  and  boundary  conditions  for  the  dynamic  shear  loading  example.  The  crack  is 
modeled  by  an  actual  discontinuity  in  the  mesh  with  a  zero  radius  crack  tip.  The  load  is  applied  as  a  velocity 
condition  that  is  ramped  up  from  0  to  16.5  m/s  in  one  microsecond  and  then  held  constant  for  the  duration  of 
the  simulation. 

The  model  parameters  are  p  =  8000  kg/m3,  E  =  190  GPa,  v  =  0.3,  Qc  =  2.213  x  104  J/m2,  k  =  0,  and 
plane  strain  is  assumed.  The  corresponding  dilatational,  shear,  and  Rayleigh  wave  speeds  are  v d  =  5654  m/s, 
vs  =  3022  m/s,  vr  =  2803  m/s.  The  length  scale  was  chosen  to  be  e  =  1.95  x  10-4  m  leading  to  a  maximum 
uniaxial  stress  of  1.07  GPa  (see  Section  2.3.1).  The  mesh  has  1024  x  1024  uniform  quadratic  elements  so 
that  h  «  e/2. 

The  simulations  were  performed  using  the  staggered  integration  scheme  described  in  Section  3.3.2  with 
the  momentum  equation  being  solved  explicitly  using  the  HHT-a  method  with  a  =  —0.1.  A  fixed  time  step 
of  At  =  1.25  x  10“  8  was  chosen,  which  is  slightly  less  than  0.9A/crit,  where  the  critical  time  step,  A/crit, 
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is  computed  as 


A  tcrit  =  (92) 

^max 

with  ccmax  the  maximum  natural  frequency  of  the  momentum  equation  determined  from  the  undamped  eigen- 
problem  and  (considering  only  the  undamped  case) 

HHT-o'  (see  Miranda,  Ferencz,  and  Hughes  (1989)) 


=  a/2(7  +  2a: (7  -  p) 
7  +  2^(7  -  p) 


(93) 


Explicit  generalized-^  (see  Hulbert  and  Chung  (1996)) 


^crit  — 


I  12(1  +  Pfr)3(2  —  pb) 

10  +  15  pb  -  pl  +  pi  -  p\ 


(94) 


The  resulting  phase-field  is  shown  in  Figure  13.  Note  that,  initially,  the  crack  starts  to  propagate  at  a  larger 
angle  then  the  angle  decreases  as  the  crack  propagates.  The  average  angle  from  the  initial  crack  tip  to  the 
point  where  the  crack  intersects  the  boundary  is  somewhat  greater  than  65°.  This  is  in  fairly  good  agreement 
with  the  experimental  results,  which  show  the  crack  propagating  at  about  70°.  The  velocity  of  the  crack  tip 
is  shown  in  Figure  15.  As  can  be  seen,  the  crack  quickly  accelerates  to  a  velocity  just  below  60%  of  the 
Rayleigh  wave  speed  and  maintains  this  velocity  until  it  decreases  as  the  crack  approaches  the  top  surface. 
Although  no  crack  tip  velocity  information  is  reported  for  the  experimental  results,  this  velocity  behavior  is 
in  good  agreement  with  behavior  reported  by  Ravi-Chandar  (1998)  for  a  similar  experiment. 

In  Figure  14  we  show  a  post-processed  plot  of  the  model  at  t  =  70  fis.  In  this  figure  we  have  scaled  the 
displacements  by  a  factor  of  5  and  removed  areas  of  the  model  from  the  plot  where  c  <  0.05  in  order  to  show 
a  representation  of  the  cracked  geometry. 


Remark:  Experience  has  shown  that  the  action  of  the  phase-field  does  not  effect  the  critical  time  step  of 
the  explicit  algorithm.  Choosing  At  <  0.9Afcr^  has  proven  sufficient  to  maintain  stability  and  is  often 
conservative,  i.e,  letting  At  =  A tcru  often  results  in  a  stable  time  step. 


4.4  Adaptive  refinement  scheme 

The  length  scale  parameter,  e,  plays  two  roles  in  the  phase-field  model:  first,  it  determines  the  width  of  the 
approximation  to  the  crack,  and  second,  as  shown  in  Section  2.3.1,  it  influences  the  magnitude  of  the  tensile 
stress  required  for  crack  nucleation.  Thus,  in  order  to  capture  fine  scale  details  of  a  crack,  or  model  materials 
with  high  nucleation  stresses,  a  small  value  for  e  is  needed.  For  example,  if  the  maximum  uniaxial  shear 
stress  is  too  low  for  the  dynamic  shear  loading  example  in  Section  4.3,  a  secondary  crack  will  nucleate  at 
the  surface  opposite  the  initial  crack  (see  Ravi-Chandar  et  al.  (2000)).  Thus,  a  small  value  of  e  is  required  to 
accurately  capture  the  crack  topology  for  this  problem.  This  in  turn  requires  a  fine  mesh  in  areas  where  the 
crack  is  located. 

To  efficiently  compute  with  fine  meshes,  as  needed  to  accurately  resolve  a  crack  for  small  values  of  e, 
we  introduce  an  adaptive  refinement  scheme.  For  this  scheme  we  choose  the  phase-field  parameter  as  a 
convenient  measure  for  determining  the  need  for  refinement.  As  has  been  shown,  the  gradients  of  the  phase- 
field  are  high  in  an  area  near  the  crack.  Away  from  the  crack  the  value  of  the  phase-field  stays  close  to  one. 
By  choosing  a  critical  threshold  of  the  phase-field  that  is  higher  than  the  value  at  which  crack  nucleation 
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(a)  t  =  20  (is 
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(b)  t  =  40  (is 


(c)  t  =  65  (is  (d)  t  =  87  (is 


Figure  13:  Evolution  of  the  crack  through  time  for  a  uniform  1024  x  1024  quadratic  NURBS  mesh  with 
1,055,242  control  points  and  e  =  1.95  x  10-4  m.  The  resulting  crack  propagation  angle  is  somewhat  greater 
than  65°  and  close  to  the  experimentally  observed  angle  of  about  70°  reported  by  Kalthoff  and  Winkler  (1987) 
and  Kalthoff  (2000). 

occurs  (c  =  0.75)  the  area  near  the  crack  is  easily  identified.  Using  a  larger  value  for  the  critical  threshold 
results  in  a  greater  area  of  refinement  (we  have  found  c  =  0.8  to  be  a  good  choice).  The  adaptive  refinement 
scheme  we  have  developed  proceeds  as  follows: 

1 .  Run  the  dynamic  simulation  to  some  termination  point 

2.  Flag  elements  where  the  phase-field  is  below  the  critical  threshold 

3.  Refine  the  flagged  elements 

4.  Rerun  the  simulation  with  the  locally  refined  mesh 
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Figure  14:  A  post-processed  plot  of  the  dynamic  shear  loading  example  at  t  =  75  /is.  The  displacements  have 
been  scaled  by  a  factor  of  5  and  areas  of  model  where  c  <  0.05  have  been  removed  from  the  plot.  Pressure  is 
measured  in  Pascals. 
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Figure  15:  Crack  tip  velocity  for  dynamic  shear  loading  example. 


5.  Repeat  steps  2 — 4  until  convergence 

4.4.1  Analysis-suitable  local  refinement  of  T-splines 

A  distinguishing  feature  of  T-splines  is  the  presence  of  T-junctions  (hanging  nodes)  in  the  mesh.  T-junctions 
maintain  locality  in  the  context  of  refinement.  This  is  in  contrast  with  NURBS  where  all  refinement  is 
global.  A  highly  localized  and  efficient  refinement  algorithm  for  analysis- suitable  T-splines  was  developed 
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Figure  16:  A  schematic  representation  of  an  adaptive  refinement  scheme  based  on  T-splines,  analysis-suitable 
local  refinement,  and  Bezier  extraction. 


by  Scott  et  al.  (2011).  This  algorithm  avoids  introducing  superfluous  control  points,  preserves  exact  geometry, 
generates  smooth  nested  spaces,  and  maintains  the  properties  of  an  analysis- suitable  space. 

In  the  context  of  isogeometric  analysis,  the  adaptive  refinement  scheme  is  based  on  Bezier  extraction  and 
analysis- suitable  local  refinement  as  shown  schematically  in  Figure  16.  First,  the  flagged  Bezier  elements  are 
used  to  determine  the  basis  functions  of  the  T-spline  that  will  be  refined.  Analysis-suitable  local  refinement 
is  then  applied  to  generate  the  refined  set  of  T-spline  basis  functions.  Bezier  extraction  is  then  applied  to  the 
refined  T-spline  to  generate  the  new  set  of  Bezier  elements. 

We  apply  this  refinement  scheme  to  the  dynamic  shear  loading  example  from  Section  4.3.  We  start  with 
a  coarse  initial  C2 -continuous  cubic  T-spline  that  has  128  x  128  Bezier  elements.  For  all  meshes,  e  is  set 
to  1.95  x  10-4m.  Elements  are  flagged  for  refinement  if  the  phase-field  parameter  is  less  than  0.8  at  any 
quadrature  point  within  the  element.  The  sequence  of  results  shown  in  Figure  17  where  each  simulation  was 
terminated  at  t  =  100  /is.  Figure  18  shows  the  sequence  of  meshes  with  the  elements  that  have  been  flagged 
for  refinement  at  the  end  of  each  iteration.  Note  that  when  the  mesh  is  too  coarse,  the  crack  propagation  is 
restricted  and  the  direction  is  incorrect.  It  is  not  until  mesh  3,  when  h  =  e,  that  the  mesh  is  fine  enough  to 
capture  the  correct  crack  path. 

Figure  19  compares  the  elastic  strain  energy  and  dissipated  energy  at  each  refinement  iteration  to  the  solu¬ 
tion  from  Section  4.3,  which  we  call  the  reference  solution.  The  elastic  strain  energy,  shown  in  Figure  19(a), 
is  over  predicted  for  the  coarse  meshes  as  a  result  of  restricted  crack  propagation.  This  plot  shows  that  it  is 
not  until  mesh  4,  when  most  of  the  element  along  the  propagation  path  are  such  that  h  =  e/2,  that  the  elastic 
strain  energy  agrees  well  with  the  reference  solution.  This  is  also  true  for  the  dissipated  energy  shown  in 
Figure  19(b). 

Table  2  lists  the  total  number  of  functions  for  each  mesh  in  the  refinement  sequence  and  the  number 
of  elements  that  were  flagged  at  the  end  of  each  simulation.  Note  that  the  final  mesh  has  53,032  cubic  basis 
functions.  This  is  compared  to  1,055,242  quadratic  basis  functions  in  the  uniformly  refined  reference  solution 
from  Section  4.3. 


Mesh  1 

Mesh  2 

Mesh  3 

Mesh  4 

Mesh  5 

Number  of  functions 
Flagged  elements 

17,755 

589 

19,992 

2,001 

27,032 

6,257 

47,824 

1,446 

53,032 

8 

Table  2:  The  number  basis  functions  before  refinement  and  the  number  of  elements  that  were  flagged  for 
refinement  for  each  mesh. 
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Figure  17:  Kalthoff  mesh  refinement  results.  Mesh  1  is  a  128  x  128  cubic  T-spline  mesh.  Bezier  elements 
were  flagged  for  refinement  if  c  <  0.8  at  any  quadrature  point  inside  the  element  and  h  =  \fa  >  1.94  x 
10-4  m  where  a  is  the  element  area. 
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Figure  18:  The  first  four  meshes  in  the  local  refinement  sequence.  The  elements  in  red  are  those  that  were 
selected  to  be  refined  at  each  step. 
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(a) 


(b) 


Figure  19:  The  (a)  elastic  strain  energy,  Jn{[(l  —  k)c 2  +  k]^  +  ^pe  }dec,  and  (b)  dissipated  energy, 
fnGc  [(c  —  l)2/ (4e)  +  e |  Vc|2]  dec,  for  the  sequence  of  refined  meshes  shown  in  Figure  17.  The  reference 
mesh  is  the  uniformly  refined  mesh  from  Figure  13. 
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4.5  Pressurized  cylinder  with  solid  elements 

A  major  benefit  of  the  phase-field  formulation  presented  here  is  that  it  extends  easily  to  three  dimensions. 
As  a  final  example,  we  show  a  three-dimensional  computation  of  a  pressurized  cylinder  with  a  spherical  end 
cap.  The  input  geometry  for  the  simulation  is  shown  in  Figure  20  where  symmetry  is  used  to  reduce  the 
computational  cost.  The  initial  crack  is  modeled  by  an  induced  phase-field  (see  Appendix  A).  A  linearly 
increasing  hydrostatic  pressure  load,  p ,  is  applied  to  the  inner  surface  as  p  =  50 1  MPa  where  t  is  the  current 
time. 


Figure  20:  Geometry  and  symmetry  conditions  for  the  pressure  vessel  simulation.  The  mesh  is  a  three- 
dimensional  thickened  T- spline. 

The  model  parameters  are  p  =  8000  kg/m3,  E  =  190  GPa,  v  =  0.3,  Qc  =  2.213  x  104  J/m2,  and 
km  0.  The  corresponding  dilatational,  shear,  and  Rayleigh  wave  speeds  are  Vd  =  5654  m/s,  vs  =  3022  m/s, 
vr  =  2803  m/s.  The  length  scale  was  chosen  to  be  e  =  2.5  x  10-3m  leading  to  a  maximum  uniaxial  stress 
of  298  MPa  (see  Section  2.3.1). 

To  construct  the  mesh,  a  C2 -continuous  cubic  T- spline  mid-surface  was  first  modeled  in  Rhino,  a  com¬ 
mercial  CAD  software  package,  using  the  T-Splines,  Inc.,  plugin.  The  initial  mid-surface  mesh  had  a  mesh 
size  of  h  «  0.01  m.  After  export,  the  surface  was  thickened  with  C1  -continuous  quadratic  functions  such 
that  there  were  eight  Bezier  elements  (eleven  functions)  through  the  thickness.  To  get  the  final  mesh,  we  used 
the  adaptive  refinement  scheme  describe  in  Section  4.4  .  The  refinement  was  applied  to  the  mid-surface  mesh 
at  each  iteration  until  h  «  e/2  in  the  area  of  the  crack.  A  new  volume  mesh  was  created  from  the  updated 
mid-surface  mesh  at  each  iteration.  The  final  mesh  is  shown  in  Figure  21.  This  mesh  contains  862,100  basis 
functions. 

The  simulations  were  performed  using  the  staggered  integration  scheme  described  in  Section  3.3.2  with 
the  momentum  equation  being  solved  explicitly  using  the  HHT-a  method  with  a  =  —0.3.  An  adaptive  time 
step  of  At  =  0.9 A Zcrit  was  used  (see  Section  4.3  for  a  definition  of  Atcrit).  To  compute  Atcrit,  we  used 
the  power  iteration  algorithm  presented  by  Benson  (1998).  The  resulting  phase-field  is  shown  at  several  time 
intervals  in  Figure  22.  A  post-processed  plot  of  the  model  at  t  =  1.76  x  10-3  s  is  shown  in  Figure  23 
with  the  displacements  scaled  by  a  factor  of  5  and  the  area  of  the  model  where  c  <  0.05  removed  from  the 
visualization. 

Figure  24  shows  two  cross-section  views  of  the  crack  at  different  times  in  the  simulation.  These  cross- 
sections  show  the  ability  of  the  phase-field  model  to  capture  three-dimensional  characteristics  of  a  crack.  We 
emphasize  that  the  computation  for  this  three-dimensional  model  did  not  require  any  additional  algorithmic 
complexity  compared  to  the  two-dimensional  models  shown  previously. 
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Figure  21:  The  final  mesh  for  the  pressurized  cylinder  example  problem.  The  volumetric  mesh  was  con¬ 
structed  by  thickening  a  mid- surface  mesh.  The  refinement  was  performed  using  the  adaptive  refinement 
scheme  describe  in  Section  4.4  which  resulted  in  a  final  mesh  containing  862,100  basis  functions. 


Remark:  This  model  demonstrates  several  key  features  of  isogeometric  analysis.  First,  the  initial  mid¬ 
surface  model  was  constructed  in  a  commercial  CAD  software  package.  This  model  was  used  directly  to 
construct  the  analysis  model,  i.e.,  there  is  no  intermediate  meshing  step.  Second,  the  smoothness  of  the  CAD 
geometry  is  represented  exactly  by  the  analysis  model.  Finally,  refinement  was  performed  directly  on  the 
CAD  model  and  the  exact  CAD  geometry  was  maintained  at  each  iteration.  This  is  illustrated  in  Figure  21. 
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(a )t=  1.02  x  10-3 s 


(b)  t  =  1.38  x  10“3  s 


(c)t  =  1.53  x  10“3s  (d)  t  =  1.76  x  10-3  s 

Figure  22:  The  results  of  the  pressurized  cylinder  example.  The  phase-field  is  shown. 
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Figure  23:  A  post-processed  plot  of  the  pressure  vessel  example  at  t  =  1.76  x  10-3  s.  The  displacements 
have  been  scaled  by  a  factor  of  5  and  areas  of  model  where  c  <  0.05  have  been  removed  from  the  plot. 
Displacement  is  measured  in  meters. 


t=  1.53xl0'3s 

Figure  24:  Cross  section  views  showing  the  three-dimensional  phase-field  profiles  of  the  crack  surfaces. 
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5  Conclusion 

We  have  extended  the  phase-field  model  for  quasi-static  brittle  fracture  presented  by  Miehe,  Hofacker,  and 
Welschinger  (2010a)  to  the  dynamic  case.  The  phase-field  model  provides  a  smooth  representation  of  a  crack 
and  removes  the  requirement  to  numerically  track  discontinuities  in  the  displacement  field.  The  width  of 
the  phase-field  approximation  of  a  discrete  crack  is  controlled  by  a  length  scale  parameter.  We  have  shown 
that  this  parameter  also  influences  the  critical  stress  at  which  crack  nucleation  occurs  and  should  therefore  be 
considered  as  a  material  parameter.  To  perform  time  integration  of  the  dynamic  model,  we  have  presented 
both  a  monolithic  and  a  staggered  time  integration  scheme.  The  staggered  scheme  provides  efficiency  and 
flexibility  in  how  the  momentum  equation  is  integrated  and  therefore  holds  greater  potential  for  large  three- 
dimensional  problems. 

We  have  studied  the  behavior  of  the  model  by  performing  numerical  experiments  for  crack  propaga¬ 
tion  and  branching.  These  examples  have  shown  that  the  phase-field  model  can  accurately  capture  complex 
dynamic  crack  propagation  behavior  in  both  two  and  three  dimensions.  In  addition,  we  have  proposed  an 
adaptive  local  refinement  strategy  that  allows  for  the  efficient  simulation  of  complex  crack  patterns.  This  re¬ 
finement  strategy  takes  advantage  of  the  character  of  the  phase-field  evolution  to  determine  when  refinement 
is  needed.  We  have  demonstrated  by  numerical  examples  the  effectiveness  of  this  refinement  scheme  in  the 
context  of  T- spline-based  isogeometric  analysis. 

As  a  finally  example,  we  have  shown  that  the  combination  of  the  phase-field  model  and  local  refinement 
strategy  provides  an  effective  method  for  simulating  fracture  in  three-dimensional  structures.  The  ability 
to  simply  and  effectively  model  three-dimensional  fracture  is  perhaps  the  most  significant  attribute  of  the 
phase-field  model. 
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A  Modeling  a  preexisting  crack  in  a  continuous  body 


For  the  numerical  examples  discussed  in  this  paper,  a  preexisting  crack  is  used  to  initialize  crack  propagation. 
The  initial  crack  is  model  as  either  a  discrete  crack  in  the  geometry  or  as  an  induced  crack  in  the  phase-field. 
For  the  induced  crack,  an  initial  strain-history  field  is  specified  such  that  an  initial  crack  in  the  phase-field  is 
defined.  To  define  the  initial  strain-history  field  we  let  l  be  a  line  that  represents  the  discrete  crack  we  wish 
to  include  and  d{x ,  l)  is  the  closest  distance  from  x  to  the  line  /.  The  strain-history  field  is  then  defined  as 


n0(x)  =  B 


d(x,l )  ^ 


d(x,  l)  <  e 
d(x,  l)  >  e 


(95) 
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where  the  magnitude  of  the  scalar  B  can  be  determined  by  letting  d  —  0  and  substituting  Ho  into  (12)2  with 
k  =  0  and  d2c/dxi  =  0  to  get 

B  =  -  —  1.  (96) 

c 

In  the  examples  presented  above,  we  have  chosen  c  =  10-3  to  be  the  value  of  the  phase-field  in  the  initial 
crack  so  that  B  =  103. 


B  Dimensionless  form  of  the  phase-field  equations 


To  improve  the  conditioning  and  scaling  of  the  fully  coupled  system  of  equations  discussed  in  Section  3.3.1, 
we  consider  the  strong  form  equations  (12)  in  their  dimensionless  form.  To  arrive  at  the  dimensionless 
formulation  we  define  a  length  scale  To  and  time  scale  To  as 


T0  = 


Gc 

ChE ’ 


To  =  T0 


(97) 


where  the  non-dimensional  constant  Ch  is  used  to  control  the  scaling  of  the  problem.  By  introducing  non- 
dimensional  space  and  time  coordinates,  and  a  non-dimensional  displacement  field  as 


x*  =  x/Lq,  t*=t/To ,  u*=u/Lo 


(98) 


we  arrive  at 


dx* 


d2i 


d(t*y 


*  \  2 


[4e*(l  -  k)H *  +  1]  c  —  4(e*) 


d2c 


(99) 


d(x*y 


=  i 


with 


a  = 


H*  = 


H 


e  = 


(100) 


ChE ’  ChE ’  T0 

In  practice,  we  have  found  that  choosing  Ch  such  that  the  size  of  the  elements  in  the  mesh  have  an  area  equal 
to  one  yields  good  results. 


Remark:  The  same  implementation  can  be  used  to  compute  with  both  the  dimensional  and  non-dimensional 
forms  of  the  equations.  If  the  dimensional  form  of  the  equations  has  been  implemented  then  the  non- 
dimensional  form  can  be  computed  by  setting  the  material  parameters  Qc  =  1,  p  =  1,  and  E  =  1  jCh\ 
scaling  the  input  geometry  by  1/Tq;  scaling  the  time  steps  by  1/Tq;  and  using  e*  in  place  of  e. 
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